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The continuous random network (CRN) model is an idealized model for perfectly coordinated amor- 
phous semiconductors. The quality of a CRN can be assessed in terms of topological and configura- 
tional properties, including coordination, bond-angle distributions and deformation energy. Using a 
variation on the sillium approach proposed 14 years ago by Wooten, Winer and Weaire, we present 
1000-atom and 4096-atom configurations with a degree of strain significantly less than the best CRN 
available at the moment and, for the first time, comparable to experimental results. 



PACS numbers: 61.43.Dq, 71.55.Jv 61.20.Ja 



I. INTRODUCTION 

The structure of amorphous semiconductors, as seen by 
theorists, is well represented by the continuous random 
network (CRN) model introduced more than 60 years 
ago by Zachariasen Q . The interest of this model lies in 
its simplicity: the only requirement of this model is that 
each atom should satisfy fully its bonding needs. In ad- 
dition, the quality of a CRN is generally determined by 
the amount of strain, as measured by the local deviations 
from crystalline environment, induced by this constraint; 
the "ideal" CRN being typically defined as that with the 
lowest spread in the bond length and bond angle distri- 
butions. 

In spite of the simplicity of the model, it has turned 
out to be difficult to actually prepare CRN realizations 
of a quality comparable to that of experiment, making it 
difficult to fully assess the real structure of amorphous 
semiconductors. The origin of this problem has gener- 
ally been attributed to weaknesses in the model-building 
community: standard approaches such as molecular dy- 
namics cannot reach time scales appropriate for full re- 
laxation. Moreover, other techniques suggest that empir- 
ical and semi-empirical potentials able to reproduce all 
properties of amorphous semiconductors are still missing 
. An alternative explanation for the fact that exper- 
imentally a lower spread in bond length and angular dis- 
tribution is observed, might be that the coordination in 
high-quality a-Si samples is significantly lower than four. 
Laaziri et al. report a coordination as low as 3.88 ||, a 
density of defects orders of magnitude higher than what 
is measured using electron-spin resonance techniques || . 
If true, this higher density of defects might easily facili- 
tate a lower spread in the bond lengths and angles, ex- 
plaining in part the discrepancy between experiment and 
theoretical models. 

Following a long tradition, one approach to shed some 



light on this discrepancy is to try to see how far it is 
possible to push the continuous random model in order 
to reach structural properties in agreement with experi- 
ment. By creating idealized networks with the same an- 
gular deviation as the experimental ones and a good over- 
all fit to the radial distribution function, it is possible to 
show that perfect coordination in amorphous silicon is 
not ruled out by the low angular deviation. This is the 
purpose of this paper which follows a long series of works 
going in the same direction [||-^| . 

Using a modified version of the Wooten- Winer- Weaire 
algorithm, we have succeeded in creating a number of to- 
tally independent 1000-atom configurations with a bond- 
angle distribution as low as 9.19 degrees, almost 2 de- 
grees below the best available numerical models with- 
out four-membered rings and on a par with experimen- 
tal values po| . The algorithm we use avoids completely 
the crystalline state, contrary to previous WWW-type 
approaches. Moreover, as shown below, the structural 
and electronic properties of the networks are excellent, 
making them ideal starting point for empirical as well as 
tight-binding or ab-initio studies jyj. 

This paper is constructed as follows. First, we review 
briefly the Wooten- Winer- Weaire algorithm and detail 
our simulation procedure. Next, we present structural 
and electronic properties of the configurations generated, 
and compare them with previous simulations and exper- 
imental results. 



II. METHODOLOGY AND DETAILS OF 
SIMULATIONS 

In the sillium approach |t]], proposed by Wooten, 
Winer and Weaire (WWW) to generate CRN structures, 
a configuration consists of the coordinates of all N atoms, 
together with a list of the bonds between them. The 
structural evolution consists of a sequence of bond trans- 
positions involving four atoms. Four atoms A, B, C, and 
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D are selected following the geometry shown in Fig. [j]; 
two bonds, AB and CD, are then broken, and atoms 
A and D reassigned, respectively, to C and B, creating 
two new bonds, AC and BD. After the transposition, all 
atoms are allowed to relax within the constraints of the 
neighbor list. 
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FIG. 1. A basic WWW event. Left: before; right: after 
the bond exchange. 

Within this approach, the generation of a CRN starts 
with a cubic diamond structure, which is then random- 
ized by a large number of bond transpositions. Af- 
ter thermalization, the network is relaxed through a se- 
quence of many more proposed bond transpositions, ac- 
cepted with the Metropolis acceptance probability 



P = Min [1, exp((£; b - E f )/k b T)] 



(1) 



where k b is the Boltzmann constant, T is the tempera- 
ture, and Ei, and Ef are the total energies of the system 
before and after the proposed bond transposition. 

The list of neighbors determines the topology, but also 
the energy of the network: independently of the distance 
between two atoms, they interact only if they are con- 
nected in the list of neighbors. With an explicit list of 
neighbors, it is possible to use a simple interaction such 
as the Keating potential p2j: 
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where a and f3 are the bond-stretching and bond-bending 
force constants, and d = 2.35 A is the Si-Si strain-free 
equilibrium bond length in the diamond structure. Usual 
values are a = 2.965eV/A 2 and (3 = 0.285 a. 

With the approach described above, along with a few 
more details that can be found in Ref. 0, Wooten and 
Weaire obtained 216-atom structures with an angular dis- 
tribution as low as 10.9 degrees. A decade later, using the 
same approach but more computing power, Djordjevic, 
Thorpe and Wooten (DTW) produced some large (4096- 
atom) networks of even better quality, with a bond angle 
distribution of 11.02 degrees for configurations without 
four-membered rings and 10.51 degrees when these rings 
where allowed 



In the present work, using a series of algorithmic im- 
provements and faster computers, we are able to generate 
structurally and electronically better networks: the 1000- 
atom configurations, for example, show a bond-angle dis- 
tribution of almost two degrees lower than DTW's model 
while our 4096-atom cell is more than one degree better. 

The improvements introduced to the sillium approach 
are the following: 

1. we start from a truly random configuration rather 
than from a molten crystalline state, thus guaran- 
teeing that the structure is not contaminated by 
some memory of the crystalline state; 

2. we evaluate the acceptance of a trial move using a 
Metropolis accept/reject procedure without doing 
full relaxation; 

3. we use a local/non-local relaxation procedure to 
limit the number of force evaluations, i.e., we re- 
lax only locally in the first ten relaxation steps af- 
ter a bond transposition (up to the third neighbor 
shell); in combination with 2), this makes the time 
per bond transposition almost independent of the 
configuration size; 

4. at regular times, we quench the structure at zero 
tem perature, with advantages outlined in section 



With these improvements, the generation of the net- 
works goes as follows. We first gene rate starting con- 
figurations as described in section II A, and quench these 



structures as described in section II C . Next, we alternate 



running at a temperature of 0.25 eV for about 100 trial 
bond transpositions per atom, and quenching. The de- 
crease in energy is almost exclusively obtained during the 
quenching, the role of the annealing at finite temperature 
is mostly to provide for a fresh starting point for the next 
quench. Once the energy is brought down to about 0.3 
eV per atom, and the angular spread around 10 degrees, 
this procedure yields diminishing returns: the annealing 
is no longer able to bring the sample to a starting point 
where the quenching leads to a lower minimum. To lower 
the energy further, we therefore also anneal the config- 
urations in different conditions for a few hundred trial 
bond transpositions per atom, like a stronger three-body 
force or a larger volume. 



A. Generating random initial CRNs 

To generate a random initial configuration, we ran- 
domly place the atoms in a box at crystalline density, 
under the constraint that no two atoms be closer than 
2.3 A. The difficult part is to connect these atoms in or- 
der to obtain a tetravalent network. We achieve this by 
starting with a loop visiting four atoms somewhere in the 
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sample, in such a way that each pair of atoms that are 
neighbors along the loop be not separated by more than 
a cut-off distance r c . This loop is gradually expanded 
until it visits each atom exactly twice; the steps of the 
loop are then the bonds in our tetravalent network. The 
expansion of the loop is achieved by randomly selecting 
a group of three atoms A, B and C, such that A is not 
four-fold coordinated and is within a distance of r c from 
B and C but not bonded to either, while B and C are 
bonded. Next, the bond BC is replaced by bonds AB 
and AC, expanding the loop by one step. This proce- 
dure is illustrated in figure ^. Initially, r c is set to some 
small value like 3 A, but then it is gradually increased 
until all atoms are four-fold coordinated. Although this 
method leads to highly strained initial configurations, it 
has the advantage that it contains absolutely no trace of 
crystallinity. 




FIG. 2. One step in the expansion of the loop, that eventu- 
ally visits all atoms four times. Three atoms A, B, and C are 
involved, of which B and C are bonded, while A is bonded to 
neither B nor C, and of which A is not four-fold coordinated. 
The bond BC is then replaced by bonds AB and AC. 

This process leads typically to CRNs whose angular 
distribution initially has a width of around thirty de- 
grees, but which reduces rapid ly to around 13 degrees 
in a single quench (see section pi C| ). In the beginning 
of this first quench, when the angular deviation is quite 
large, sometimes a pair of atoms is closeby without be- 
ing bonded; to eliminate such artefacts that result from 
the fact that within the Keating potential atoms only in- 
teract if they are explicitly bonded, we replace a bond 
of each of these atoms by a bond between these atoms 
and another bond between their neighbors (conserving 
four- fold coordination) . 

B. Avoiding complete relaxation of trial 
configurations 

In the standard sillium approach, a move consists of a 
bond transposition followed by full structural relaxation 
and an accept /reject step according to the Metropolis 
criterion ([!]). An alternative implementation is to first 
decide a threshold energy given by 

E t = E b -k b Tlog(s), (3) 



where s is a random number between and 1. The pro- 
posed move is then accepted only if Ef < E t . This pro- 
cedure is exactly equivalent to the usual Metropolis pro- 
cedure. By fixing the threshold for Ef before a trans- 
position is attempted, it is however possible to reject 
the move as soon as it becomes clear that this thresh- 
old cannot be reached, i.e., before the configuration is 
fully relaxed. Since the energy is harmonic around the 
minimum, the decrease in energy obtained during fur- 
ther relaxation is approximately equal to the square of 
the force times some proportionality constant c/, so that 
during the relaxation the final energy can be estimated 
to be 

Ef^E-c f \F\\ (4) 

If, at any moment during the relaxation, E—Cf\F\ 2 > E t , 
the trial move is rejected and a new one is started. Such 
a method requires, of course, a conservative choice for 
Cf, in our units, the proportionality constant Cf in well- 
relaxed configurations is always well below 1. To account 
for anharmonicities, we do not reject any move during the 
first five steps of relaxation. 

Since much less than one percent of the proposed 
moves are accepted in well-relaxed samples, avoiding 
spending time on moves that are eventually rejected can 
produce a significant gain in efficiency; using this im- 
provement, we observed a speed-up of close to an order 
of magnitude, so that all bond transpositions in a 1000- 
atom network can be attempted in less then 3 minutes 
on a 500 Mhz DEC- Alpha workstation. 

C. Efficient quenching 

Further optimizations are possible in the case of zero 
temperature. Since the threshold energy (Eq. (§)) is 
constant, a proposed bond transposition that is once re- 
jected will keep on being rejected, as long as no other 
bond transpositions are accepted in the mean time. A 
combination of four atoms ABCD with bonds AB, BC 
and CD can be selected in N ■ 4 • 3 • 3/2 times, so there 
are 187V possible bond transpositions. We mark all bond 
transpositions that were rejected since the last accepted 
bond transposition to avoid retrying these. Once all bond 
transpositions have been tried but rejected, the quench- 
ing is complete. At this stage, the system is not only 
at a local energy minimum (i.e. a point in phase space 
where the force is zero and all eigenvalues of the hessian 
are positive), but no single bond transposition can lower 
the energy. The configurations we discuss here, have this 
property. 

In the standard sillium approach, the creation of four- 
membered rings is disallowed. Following DTW we 
find that especially for quenching the relaxation is sig- 
nificantly helped by allowing for four-membered rings, 
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because of the large extra number of pathways accessible 
to the system. At the end of the quenching, the few four- 
membered rings that are created can easily be removed 
one by one, by chosing the energetically most favorable 
bond transposition in which bond AB is part of the four- 
membered ring (and where no new four-membered rings 
are introduced). Typically, the energy increases by less 
than an eV per removed four-membered ring. 

III. RESULTING CONFIGURATIONS 

We present here results for three different configura- 
tions: two 1000-atom cells and one with 4096 atoms. In 
Table ||, we compare our configurations, relaxed with the 
Keating potential used in the modified WWW algorithm, 
with those of Djordjevic, Thorpe and Wooten j|. We also 
provide the irreducible ring statistics, i.e., those rings in 
which no two atoms are connected through a sequence of 
bonds that is shorter than the sequence along the ring. 
We also provide the ring statistics for all n-membered 
rings in order to compare with Ref. M and other pa- 
pers in the literature (Table |l|) . Table [j shows that the 
strain per atom in our structures is significantly below 
that of DTW. One of the standard measurements to eval- 
uate the quality of a model is the coordination number 
as computed based on the radial distribution function 
(RDF). Using the minimum of the RDF between the first- 
and second-neighbor peak and after relaxation with the 
Keating model, the first two configurations are perfectly 
tetravalent. The 4096-atom configuration has 0.1% of 
5-fold defects. 

TABLE I. Energetic and structural properties of models 
relaxed with the Keating potential. The first two models, 
DTW (1) and DTW (2) are the models prepared in [§] and re- 
fer, respectively, to a model with and without four-membered 
rings. Configurations 2 and 3 are 1000-atom models prepared 
according to the procedure described here and '4096' is a 
4096-atom model prepared the same way. All three models 
are without four-membered rings. The ring statistics are for 
irreducible rings and po is based on ro = 2.35 A. 





DTW (1) 


DTW (2) 


Conf. 2 


Conf. 3 


4096 


E(eV) /atom 


0.336 


0.367 


0.267 


0.264 


0.304 


p/Po 


1.000 


1.000 


1.043 


1.040 


1.051 


{r)/r 


0.996 


0.997 


0.982 


0.982 


0.980 


Ar/ro (%) 


2.52 


2.65 


3.94 


0.371 


4.17 


(9) 


109.24 


109.25 


109.30 


109.27 


109.28 


A8 


10.51 


11.02 


9.21 


9.20 


9.89 


Rings/atom 












4 


0.015 


0.000 


0.000 


0.000 


0.000 


5 


0.491 


0.523 


0.472 


0.480 


0.490 


6 


0.698 


0.676 


0.761 


0.750 


0.739 


7 


0.484 


0.462 


0.507 


0.515 


0.467 


8 


0.156 


0.164 


0.125 


0.116 


0.148 


9 






0.034 


0.033 


0.035 



Another important quantity that can be compared 
with experiment is the width of the bond angle distri- 
bution. Experimentally this quantity can be extracted 
from the radial distribution function, or from the Ra- 
man spectrum — using a relation proposed by Beeman 
et al. Jl3[ |. The most recent measurement, taken on an- 
nealed samples prepared by ion bombardment and us- 
ing the second-neighbor peak of the radial distribution 
function, gives 10.45 and 9.63 degrees, respectively, for 
as-implanted and annealed samples [ fl0| ]. Our configura- 
tions, relaxed with the Keating potential, present there- 
fore a bond angle distribution slightly narrower than ex- 
perimental values. (This is to be expected of the "right" 
structure since the theoretical models are taken at zero 
K.) 

While structural averages provide a good idea of the 
overall quality of a model, they do not say much regard- 
ing the local environments. It is therefore important to 
look also at the electronic properties of these models: 
even small densities of highly strained geometries or de- 
fect atoms will be picked up as states in the gap of the 
electronic density of states (EDOS). In the last few years, 
it has become possible to compute the electronic struc- 
ture of multi-thousand-atom configurations. Here, we 
show the electronic density of state for Configuration 2, 
a 1000-atom configuration. Because of the costs of do- 
ing a full ab-initio atomic relaxation, we have relaxed 
the cell with the Keating potential and used the Fireball 
local-basis ab-initio code to obtain the electronic density 
of state [^4). Previous work showed that configurations 
relaxed with a Keating potential demonstrated little fur- 
ther relaxation with the Fireball code |L5| so that the 
results presented here are unlikely to change very much 
during further relaxation. Figure I shows the EDOS 
smoothed with a gaussian of width 0.01 eV. A remark- 
able feature of this EDOS is the absence of states in the 
gap, leading to a perfect gap of 1.3 eV. The generation of 
defect-less models is very important for our understand- 
ing of the electronic dynamics and the role of defects in 
disordered semiconductors. The decay of the valence tail, 
the Urbach tails, can be reasonably well approximated by 
an exponential — p(E) oc exp(— E/Eq) — with Eq = 0.2 
eV, in agreement with previous calculations |jl5|| . 

Although we get good structures using the Keating 
potential, it is important to verify the stability of these 
networks when relaxed with a more realistic interaction 
potential that does not require a pre-set list of neighbors. 
There exists no empirical potential at the moment that 
can describe fully the properties of a-Si; we use a modified 
Stillinger- Weber potential, where the three-body contri- 
bution to the energy is enhanced by 50 % with respect to 
the two-body term. This ad-hoc modification was shown 
to produce good structural properties for amorphous sil- 
icon 

After relaxation at zero pressure, the two 1000-atom 
configurations remain perfectly coordinated. The 4096- 
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atom cell, less well relaxed, develops a few coordination 
defects : based on the first minimum in the RDF, it 
presents respectively 0.4% and 0.3% of three-fold and 
five-fold coordinated atoms. 

Table [0| presents the structural and energetic proper- 
ties of the relaxed CRNs at zero pressure. For all con- 
figurations, the bond angle distribution widens and the 
density decreases significantly compared to the Keating- 
relaxed structures. For the 1000-atom configurations, the 
local relaxation with the modified Stillinger- Weber po- 
tential did not result in a change of topology and the 
total energies are very low compared with previous mod- 
els [[[7). We therefore expect that the structures will be 
stable with any reasonable potential. 
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FIG. 3. Top: Electronic density of states for the 
1000-atom model 2 as obtained using ab-initio tight-binding. 
Bottom: Close-up on the gap region. The dashed curve is an 
exponential fit, p(E) cx exp(—E/Eo), with Eo = 0.2eV. 

Figure |] shows a comparison of Configuration 3 with 
experimental data obtained by Laaziri et al. on annealed 



a- Si samples prepared by ion bombardment. The agree- 
ment between the two is excellent except for some dis- 
crepancy in the height of the third-neighbor peak. Such 
an agreement must only be seen as a sign that the topol- 
ogy might be right, however: configurations differing 
widely in their topology can easily produce similar radial 
distribution functions |Q . The same figure also presents 
the bond-angle distribution for Configuration 3 relaxed 
with both the Keating and the modified Stillinger- Weber 
potentials. As expected for a perfectly coordinated sam- 
ple, the distribution is smooth and presents a single peak 
centered at the tetrahedral angle. 
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FIG. 4. Top: Radial distribution function for Configura- 
tion 3 relaxed with the modified Stillinger- Weber potential. 
Solid line: experimental results from Ref. [[u)|. To obtain a 
better fit, the computer-generated structure is scaled by a 
linear factor of 0.99. Bottom: Bond angle distribution for 
Configuration 3 relaxed with Keating and the modified Still- 
inger- Weber potential. The curve is smoothed with a gaussian 
of width 2 degrees. 
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To compare with previous molecular dynamics stud- 
ies |H| , we have also relaxed our cells with the standard 
Stillinger- Weber potential, which is known to give an in- 
correct amorphous structure. After relaxation of Config- 
uration 2, we find 17 atoms with five-fold coordination 
and three three- fold coordinated ones; similar results are 
found with the two other networks. The resulting config- 
urational energy, given in Table || compares favorably 
with molecular dynamical results. 

Contrary to the topological properties, which seem rel- 
atively independent of the details of the potential, we 
see that the ideal density of amorphous silicon compared 
with that of the crystal changes qualitatively as a func- 
tion of the potential used. Configurations relaxed at zero 
pressure with the Keating potential show a reproducible 
dcnsification by 2 % while the modified Stillinger- Weber 
potential, also at zero pressure, leads to a structure which 
is up to 6 % less dense. The latter results are in qualita- 
tive agreement with experiment [pp[ and previous simu- 
lations using a similarly modified potential f| . Recently, 
Laaziri and collaborators have pointed to the lower den- 
sity of a-Si as an explanation for the relatively low coor- 
dination measured by x-ray scattering. Our results, on 
the contrary indicate that there is very little dependence 
between the density of the amorphous material and its 
topology, at least within the application of our two em- 
pirical potentials. A volume change at the percent level 
should therefore have very little impact on the topology 
and will reflect mostly some fine details of the real atomic 
interactions. 

TABLE II. Structural properties of our configurations after 
relaxation with the modified Stillinger- Weber (mSW) poten- 
tial. Except for the 4096-atom configuration, the topology has 
remained unchanged (see text.) For comparison purposes, the 
total ring number per atom (including reducible ones) is also 
reported, as well as the energy after relaxation with the orig- 
inal Stillinger- Weber (SW) potential; MD-prepared samples 
give -4.088 eV/atom M. 





Sample 2 


Sample 3 


4096 




E(eV)/atom (mSW) 


-4.026 


-4.034 


-3.990 




E(eV)/atom (SW) 


-4.126 


-4.133 


-4.106 




p/Po 


0.947 


0.950 


0.936 




{r)/r 


1.018 


1.017 


1.020 




Ar/ro (%) 


2.9 


2.7 


0.032 






109.25 


109.24 


109.20 




M 


9.77 


9.70 


10.51 




Rings/atom 








4 


0.000 


0.000 


0.001 




5 


0.472 


0.480 


0.489 




6 


0.840 


0.847 


0.830 




7 


1.011 


1.023 


0.979 




8 


2.025 


2.002 


2.064 





IV. CONCLUSIONS 

We have presented here modifications to the Wooten- 
Winer-Weaire algorithm that allows one to produce large 
high-quality continuous random networks without pass- 
ing at all by a crystalline phase. Structural and elec- 
tronic properties of the networks produced are excellent 
and they compare favorably with experiment. 

The coordinates of the three configurations discussed 
here, as well as a 10 000-atom sample under preparation, 
are available by request. 
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